!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! This module collects all the utilities connected with the
!! definition of the stabilization parameter \f$\tau\f$ for the LDG-H
!! method.
!!
!! \n
!!
!! In general, it is allowed to specify a stabilization parameter
!! regrdless of the type of the FE basis. Use <tt>stab_method =
!! 'none'</tt> to avoid stabilization.
!! 
!! The \f$\tau\f$ parameters are collected in the vector \c taus,
!! which follows the order of the elements in the grid.
!! 
!! For the sake of generality, the parameter \f$\tau\f$ is allowed to
!! vary along each side, and is computed at the boundary quadrature
!! nodes. For this reason, this module also depends on \c mod_base and
!! the chosen basis.
!<----------------------------------------------------------------------
module mod_tau

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave

 use mod_output_control, only: &
   mod_output_control_initialized, &
   base_name

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me

 use mod_grid, only: &
   mod_grid_initialized, &
   t_v, t_s, t_e,       &
   p_t_v, p_t_s, p_t_e, &
   t_grid,              &
   affmap,              &
   diameter

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_testcases, only: &
   mod_testcases_initialized, &
   coeff_diff,  &
   coeff_adv

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_tau_constructor, &
   mod_tau_destructor,  &
   mod_tau_initialized, &
   t_eltau,  &
   taus,     &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 
 !> \f$\tau\f$ parameter for one element
 type t_eltau
   integer :: s_tau !< local index of side \f$s_\tau\f$
   !> local indexes of non-tau sides: \f$\partial K\backslash s_\tau\f$ (d)
   integer, allocatable :: s_ntau(:)
   !> values \f$\tau\f$ at nodes \c xigb (ms,d+1)
   real(wp), allocatable :: tau(:,:)
 end type t_eltau

! Module variables

 ! public members
 type(t_eltau), allocatable :: taus(:) !< element parameters (ne)

 logical, protected ::               &
   mod_tau_initialized = .false.

 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_tau'

 interface write_octave
   module procedure write_tau_struct
 end interface

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_tau_constructor(stab_method,tau_scaling,bartau_scaling, &
                                grid,base,write_grid)
  character(len=*), intent(in) :: stab_method
  integer, intent(in) :: tau_scaling, bartau_scaling
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base
  logical, intent(in) :: write_grid

  integer :: ie, is, l, idx(grid%d+1), fu, ierr
  real(wp) :: h, tau_side(grid%d+1),                         &
    xgb(grid%m,base%ms), mu(grid%m,grid%m,base%ms,grid%d+1), &
    adv(grid%m,base%ms,grid%d+1), diff_ave, n(grid%m), adv_ave
  character(len=100+len(stab_method)) :: message
  character(len=*), parameter :: out_file_suff = '-tau.octxt'
  character(len=len(trim(base_name))+len(out_file_suff)) :: out_file_name
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_fu_manager_initialized.eqv..false.).or.&
       (mod_octave_io_initialized.eqv..false.).or. &
       (mod_master_el_initialized.eqv..false.).or. &
       (mod_grid_initialized.eqv..false.)     .or. &
       (mod_base_initialized.eqv..false.)     .or. &
       (mod_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_tau_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   idx = (/(is, is=1,grid%d+1)/)
   allocate(taus(grid%ne))

   iedo: do ie=1,grid%ne
     ! Choose a first s_tau randomly. This index will be
     ! double-checked a posteriori, to ensure that tau is maximum on
     ! s_tau.
     taus(ie)%s_tau = select_random_stau(ie,grid%d)
     allocate(taus(ie)%s_ntau(grid%d))
     taus(ie)%s_ntau = pack(idx,idx.ne.taus(ie)%s_tau)
     allocate(taus(ie)%tau(base%ms,grid%d+1))

     stab: select case(trim(stab_method))
      case('none')

       ! tau is set to zero
       taus(ie)%tau = 0.0_wp

      case('upwind')
       
       ! 0) evaluate problem coefficients: diffusion and advection at
       ! the boundary Gaussian nodes
       do is=1,grid%d+1 ! side loop
         xgb = affmap(grid%e(ie),base%xigb(:,:,is))
         mu(:,:,:,is) = coeff_diff(xgb)
         adv(:,:,is)  = coeff_adv(xgb)
       enddo

       ! 1) compute the elliptic tau
       ! we need the mean diffusion coefficient on e_tau
       diff_ave = 0.0_wp
       ! get the normal n (the orientation does not matter)
       n = grid%e(ie)%n(:,taus(ie)%s_tau)
       do l=1,base%ms
         diff_ave = diff_ave + base%wgb(l,taus(ie)%s_tau)* &
           dot_product(n,matmul(mu(:,:,l,taus(ie)%s_tau),n))
       enddo
       diff_ave = diff_ave/sum(base%wgb(:,taus(ie)%s_tau))
       h = diameter(grid%e(ie)%s(taus(ie)%s_tau)%p)
       taus(ie)%tau(:,taus(ie)%s_tau) = diff_ave/h

       ! 2) add the hyperbolic contribution
       taus(ie)%tau(:,taus(ie)%s_ntau) = 0.0_wp
       do is=1,grid%d+1
         ! we need the mean flow velocity 
         adv_ave = 0.0_wp
         ! outward normal
         n = grid%e(ie)%n(:,is)
         do l=1,base%ms
           adv_ave = adv_ave + base%wgb(l,is)* &
             dot_product(adv(:,l,is),n)
         enddo
         adv_ave = adv_ave/sum(base%wgb(:,is))
         if(adv_ave.lt.0.0_wp) & ! inflow
           taus(ie)%tau(:,is) = taus(ie)%tau(:,is) - adv_ave
       enddo

      case('scaling')

       ! define tauK
       h = diameter( grid%e(ie)%s(taus(ie)%s_tau)%p )
       taus(ie)%tau(:,taus(ie)%s_tau) = h**tau_scaling

       ! define bar_tauK
       bartauk: if(bartau_scaling.gt.10) then ! set bar_tauK to zero
         do is=1,grid%d
           taus(ie)%tau(:,taus(ie)%s_ntau(is)) = 0.0_wp
         enddo
       else
         do is=1,grid%d
           h = diameter( grid%e(ie)%s(taus(ie)%s_ntau(is))%p )
           taus(ie)%tau(:,taus(ie)%s_ntau(is)) = h**bartau_scaling
         enddo
       endif bartauk

      case default
        write(message,'(A,A,A)') &
          'Unknown stabilization method "',trim(stab_method),'".'
        call error(this_sub_name,this_mod_name,message)
     end select stab

     ! Check that indeed tau is maximum on s_tau, and if necessary
     ! redefine s_tau
     tau_side = sum(taus(ie)%tau,1)/real(base%ms,wp)
     if(tau_side(taus(ie)%s_tau).lt.maxval(tau_side(taus(ie)%s_ntau))) then
       ! we need to redefine s_tau
       taus(ie)%s_tau = maxloc(tau_side,1) ! use dim to return a scalar
       taus(ie)%s_ntau = pack(idx,idx.ne.taus(ie)%s_tau)
     endif

   enddo iedo

   ! write the octave output
   if(write_grid) then
     out_file_name = trim(base_name)//out_file_suff
     call new_file_unit(fu,ierr)
     open(fu,file=trim(out_file_name), &
          status='replace',action='write',form='formatted',iostat=ierr)
       call write_octave(taus,'taus',fu)
     close(fu,iostat=ierr)
   endif

   mod_tau_initialized = .true.

 end subroutine mod_tau_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_tau_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_tau_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(taus)

   mod_tau_initialized = .false.
 end subroutine mod_tau_destructor

!-----------------------------------------------------------------------
 
 !> Randomly select \f$s_\tau\f$. We don't use the standard random
 !! number generator to ensure that the same side is chosen with
 !! different compilers.
 !<
 pure function select_random_stau(ie,d) result(s_tau)
  integer, intent(in) :: ie, d
  integer :: s_tau
 
  real(wp), parameter :: pi = 3.14159265358979_wp
  real(wp) :: nu, staun

   nu = 1.0_wp/(2.0_wp*pi)*91.0_wp/50.0_wp
   staun = sin(2.0_wp*pi*nu*real(ie,wp))
   
   ! floor to {1,2,...,d+1}
   staun = (real(d+1,wp)-1e-6_wp)*((1.0_wp/2.0_wp)*staun+0.5_wp)+1.0_wp
   
   s_tau = floor(staun);
 
 end function select_random_stau
 
!-----------------------------------------------------------------------
 
 subroutine write_tau_struct(tau,var_name,fu)
  integer, intent(in) :: fu
  type(t_eltau), intent(in) :: tau(:)
  character(len=*), intent(in) :: var_name
 
  integer :: i
  character(len=*), parameter :: &
    this_sub_name = 'write_tau_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 3' ! number of fields

   ! field s_tau
   write(fu,'(a)')      '# name: s_tau'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(tau)
   do i=1,size(tau)
     call write_octave(tau(i)%s_tau,'<cell-element>',fu)
   enddo
   ! field s_ntau
   write(fu,'(a)')      '# name: s_ntau'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(tau)
   do i=1,size(tau)
     call write_octave(tau(i)%s_ntau,'r','<cell-element>',fu)
   enddo
   ! field tau
   write(fu,'(a)')      '# name: tau'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(tau)
   do i=1,size(tau)
     call write_octave(tau(i)%tau,'<cell-element>',fu)
   enddo

 end subroutine write_tau_struct
 
!-----------------------------------------------------------------------

end module mod_tau

